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Abstract 

Paying attention to the difference of density of states, Aln g(E) = In g(E + AE) — In g(E), 
we study the convergence of the Wang-Landau method. We show that this quantity is a good 
estimator to discuss the errors of convergence, and refer to the 1/t algorithm. We also examine 
the behavior of the lst-order transition with this difference of density of states in connection with 
Maxwell's equal area rule. A general procedure to judge the order of transition is given. 
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The Monte Carlo simulation has become a standard method to study many-body prob- 
lems in physics. However, we sometimes suffer from the problem of slow dynamics in the 
original Metropolis algorithm [1]. One attempt to conquer the problem of slow dynamics is 
the extended ensemble method; one uses an ensemble different from the ordinary canonical 
ensemble with a fixed temperature. The multicanonical method [2, 3], the parallel tem- 
pering, or the exchange Monte Carlo method [4, 5] and the Wang-Landau (WL) algorithm 
[6] are examples. The WL method is an efficient algorithm to calculate the energy density 
of states (DOS), g(E), with high accuracy, and was successfully applied to many problems 
[7, 8]. The refinement and convergence of the WL method were argued [9, 10], but the 
convergence property is still a topic of discussions [11]. The search for optimal modification 
factor was discussed [12], and in connection with the WL method, the 1/t algorithm [13, 14] 
was proposed. Moreover, tomographic entropic sampling scheme has been proposed as an 
algorithm to calculate DOS [15]. 

In this paper, we investigate the convergence properties of the WL method, paying special 
attention to the difference of DOS. We argue its relevance to the lst-order transition. We 
provide a general strategy to judge the order of transition. 

Let us briefly review the WL algorithm. A random walk in energy space is performed 
with a probability proportional to the reciprocal of the DOS, l/g(E), which results in a 
flat histogram of energy distribution. Actually, we make a move based on the transition 
probability from energy level E\ to E 2 

p( El ^E 2 ) = mm\l,^}. (1) 
L g(E 2 )i 

Since the exact form of g{E) is not known a priori, we determine g{E) iteratively. Intro- 
ducing the modification factor fa, g(E) is modified by 

In In + In £ (2) 

every time the state is visited. At the same time the energy histogram h(E) is updated as 

h(E) ->■ h(E) + 1. (3) 

The modification factor /j is gradually reduced to unity by checking the 'flatness' of the 
energy histogram. The 'flatness' is checked such that the histogram for all possible E is not 
less than some value of the average histogram, say, 80%. Then, /$ is modified as 

In f i+1 = l - In f % (4) 
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and the histogram h(E) is reset. As an initial value of fi, we choose fo = e; as a final value, 
we choose In/, = 2~ 26 , that is, / 2 6 — 1.000 000 01, for example. 
We first treat the Ising model, whose Hamiltonian is given by 



Here, J is the coupling and Oi is the Ising spin (±1) on the lattice site i. The summation 
is taken over the nearest neighbor pairs Periodic boundary conditions are employed. 

Throughout this paper, we measure the energy in units of J unless specified; in other words, 
we put J — 1. 

We calculate \ng(E) with the use of the WL method, and consider the difference of 
In g(E), which is defined as 



For the Ising model, AE = 4 J. The exact value of g{E) for the two-dimensional (2D) Ising 
model is available due to Beale [16]. The deviation of the calculated value of A\n g(E) from 
the exact value of Beale [16] can be used as a measure of the accuracy of the calculation. 

We plot the overall behavior of A\n g(E) for the 2D Ising model with system size L = 
32 in Fig. 1. The data for the modification step % = 14, 18 and 22 are given for a single 
measurement. In the accuracy of this plot, little difference in % is appreciable except for 
small and large E. The enlarged plot near E = is given in the inset of Fig. 1, and the 
data for % = 14, 18 and 22 are compared to the exact value of Beale [16]. We see that 
the calculated value of A In g(E) approaches the exact value as the modification factor fi 
approaches 1. The deviation becomes smaller as % increases. The advantage of using Eq. (6) 
is that we can directly discuss the error of DOS without caring about the normalization of 
g{E). Since the transition probability depends on the difference of Ing(Ei) and In g(E2), 
this quantity of difference is essential in the method calculating the energy DOS compared 
to g(E) itself. We note that the quantity of difference was also used in the argument of 
accuracy and convergence of the WL method by Morozov and Lin [17]. 

To see the convergence of errors more explicitly, we consider the total sum of the squared 
error of A\n g(E) - A \ng(E) exact ; 




(5) 



A\ng(E) = In g(E + AE) — In g(E). 



(6) 
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2JN-12J 



A 2 = 



N-A 



(Alns(£)-Alns(£) cxact ) 2 



(7) 



E=-2JN+8J 
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FIG. 1: (Color online) Plot of Alng(E) for the 2D Ising model with L = 32. Data for i = 14, 18 
and 22 are given. In the inset, the enlarged plot near E = is shown. The exact value due to 
Beale [16] is also given in the inset for comparison. 

For the 2D Ising model, we note that A In g(-2JN + 8 J) cxact = - A In g(2JN - 12 J) cxact = 
In 2. 

In Fig. 2, we plot A 2 , Eq. (7), as a function of the modification step i up to 26 for L = 
32. The average is taken for 10 samples. We see that A 2 becomes smaller with the increase 
of i. However, the errors are saturated even though we repeat the iteration process up to 
% — 26. Such saturation of convergence of the WL method was pointed out by Yan and de 
Pablo [18]. To overcome this difficulty, a modified version of the WL algorithm in which 
the refinement parameter is scaled down as 1/t (with t the Monte Carlo time) was proposed 
[13, 14]. It is interesting to compare the performance of the 1/t algorithm and that of the 
original WL method in this quantity of difference of DOS. In the 1/t algorithm, starting from 
the same condition as the original WL algorithm, the modification factor In /j is reduced as 
1/t instead of checking the flatness condition after the condition ln/j < 1/t is satisfied. The 
final value of In / should be fixed from the outset. In Fig. 2, we also plot the data for the 
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FIG. 2: (Color online) Convergence of errors, A 2 , for the 2D Ising model with L = 32 as a function 
of the modification step i. The convergence of the original WL algorithm is compared with that 
of 1/t algorithm. In the 1/t algorithm after the rule of modification is changed, the meaning of i 
is such that MCS is 2\ 

1/t algorithm. In the case of L = 32, the modification process is changed from the original 
WL scheme to the 1/t one around i =21 or 22. In the range of 1/t scheme the actual MCS 
is fixed as 2*, which is different from the case of the original WL scheme. We clearly confirm 
the efficiency of the 1/t algorithm. In the discussion of the convergence of 1/t algorithm, 
the quantity \ng(E) — \ng(E gronn a) was used [11]. The quantity given by Eq. (6) is more 
flexible as it can be treated even if the ground state of the system is unknown as spinglass 
problems. 

We can consider the deviation from the exact value, as in Eq. (7), for the 2D Ising model. 
In order to investigate the convergence behavior of the system whose exact g(E) is not 
available, we may employ another strategy. For example, we may consider the relative error 
of the data for i and those for i — We leave the detailed analysis to a separate publication. 

Next we deal with the 2D ten-state Potts model, which is a typical model to exhibit the 
lst-order transition. This model was used to show the effectiveness of the multicanonical 
simulation by correctly estimating the interfacial free energy [3], which was later proved by 



5 



0.5 



1 



E/N 



FIG. 3: (Color online) Plot of Aln g(E) of the 2D ten-state Potts model for L = 64 (upper) and 
L = 128 (lower) as a function of E/N . The data for the modification factor /j with i = 14, 18 and 
22 are given. 

the explicit formula [19]. The Hamiltonian of the g-state Potts model is given by 



Here, Si is the Potts spin which takes 1, • • ■ , q. We note that for q = 2 the Potts model 
becomes the Ising model, although the unit of J in Eq. (8) for the Potts model is twice as 
J in Eq. (5) for the Ising model. 

We plot the difference of DOS, Eq. (6), of the 2D ten-state Potts model in Fig. 3. The 
data for L = 64 (upper) and those for L = 128 (lower) are given as a function of E/N. We 
show how the data converge as i increases by giving the data for i = 14, 18 and 22 with a 
single measurement. We clearly see the convergence of errors with the increase of i. 




(8) 
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FIG. 4: (Color online) Schematic illustration of Maxwell's equal area rule. 



The systems which show the lst-order transition have double maximum structure in the 
thermodynamic limit at the lst-order transition temperature T c when we plot the free energy 
— /3F = lng(E) — (3E as a function of E. Then, Alng(E), which is defined as Eq. (6), has 
an S-like structure with minimum and maximum. We clearly find this structure in Fig. 3. 
We note that the overall size dependence is small in this plot, but the detailed analysis is 
given later. 

The lst-order transition temperature, T c = 1/(3 C , can be estimated by Maxwell's rule as 
in thermodynamics. A schematic illustration of Maxwell's rule is shown in Fig. 4. The value 
of /3, which separates the shaded region and gives the same area, becomes the lst-order 
transition temperature /3 C . This equal area rule is proved by the following. The condition 
that the two areas of the shaded region are equal is given by 



which leads to the condition that the double maxima in In g(E) — f3E take the same value. 
In the thermodynamic limit, the difference A\n g(E) becomes the differential d\ng{E)/dE. 
The area of the shaded region, Eq. (9), is related to the interfacial free energy [3, 19]. 

To see the S'-like structure explicitly, we make an enlarged plot along y-axis of A In g(E) 
for L = 64 (upper) and 128 (lower) in Fig. 5. The modification step i is 22. In this plot we use 
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FIG. 5: (Color online) Enlarged plot of A \ng(E) of the 2D ten-state Potts model for L = 64 (upper) 
and 128 (lower). The modification step i is 22. The smoothed values with moving-average method 
are given. The lst-order transition temperature /3 C = ln(l + \/To) = 1.42606 in the thermodynamic 
limit is also shown by straight line for convenience. 

the data with the smoothing process, (f(E-2AE)+4*f(E-AE)+6*f(E)+4:*f(E+AE) + 
f{E + 2A£))/16 with f(E) = A In g(E), to reduce fluctuations. For the 2D ten-state Potts 
model, the lst-order transition temperature is given by j3 c = ln(l + Vl~0) = 1.42606. We 
give this value in Fig. 5 for convenience; we see that Maxwell's rule works. We can estimate 
P c and the interfacial free energy from the S'-like curve for each size. We observe the size 
dependence in Fig. 5; the area of the shaded region illustrated in Fig. 4 is proportional to 
1/L, which reflects on the finite size scaling of the lst-order transition. 

We may provide a general strategy to judge the order of transition for any system. We 
plot A In g(E) and check whether there is an S'-like structure. If the system shows the 1st- 
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order transition, we can locate the transition temperature by Maxwell's rule. The behavior 
of the lst-order transition can be observed in the early stage of WL iteration, that is, for 
small %. If we investigate In g(E) — (3E as in usual way, we have to search for /3 which gives 
the same value for two maxima. 

To summarize, we have shown that the difference of h\g(E) is a good quantity for the 
WL method. Less attention has been given to the quantity A In g(E) so far, although some 
efforts were made in the discussion of accuracy and convergence of the WL method [17]. 
Comparing with the exact value of the 2D Ising model, we have shown the convergence 
property of the WL method. That is, we have shown how errors become smaller for larger 
i, where i is the step of the modification factor fi for the criterion of 'flatness' condition. 
We have confirmed the efficiency of the 1/t algorithm; we have shown that the quantity 
A In g(E) is a good estimator for the analysis of errors of the simulation method to calculate 
the energy DOS. 

We have also shown that A In g{E) is a good estimator for the lst-order transition. We 
have investigated the 2D ten-state Potts model. The lst-order transition is observed in the 
S'-like behavior of Aln g(E). We have shown that Maxwell's equal area rule determines 
the lst-order transition temperature. Although the statement is rigorously realized in the 
thermodynamic limit, we observe the behavior of the lst-order transition even for small 
system size and for small i of the modification step. We assert that we provide a general 
procedure to study the order of transition for any system. 

The extension of this calculation to continuous spin models is straightforward [20]. The 
application to quantum Monte Carlo simulation [21] for checking the order of transition is 
highly desirable. The application to first-principle calculation of electric structure [22] and 
to protein systems [23] may be other interesting topics. s 

Before closing, we mention about the calculation techniques. We have used the parallel 
calculation with multiple random walkers for the WL algorithm using the GPU (graphic 
processing unit) with CUDA (common unified device architecture). The details of the GPU- 
based calculation will be given elsewhere. 

We thank Tasrief Surungan for valuable discussions. This work was supported by a 
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